Dynamics of a flexible polymer in planar mixed flow 

Dipjyoti Das^, Sanjib Sabhapandit^ and Dibyendu Das^ 

^ Department of Physics, Indian Institute of Technology Bombay, Powai, Mumbai - 400076, 
India 

^ Raman Research Institute, Bangalore 560080, India 

E-mail: dipjyotiSphy . iitb . ac . in, sanjibOrri .res . in, dibyenduSphy . iitb . ac . in 

Abstract. We present exact spatio-temporal correlation functions of a Rouse polymer chain 
submerged in a fluid having planar mixed flow, in the steady state. Using these correlators, 
determination of the time scale distribution functions associated with the first-passage tumbling 
events is difficult in general; it was done recently in Phys. Rev. Lett. 101, 188301 (2008), for 
the special case of "simple shear" flow. We show here that the method used in latter paper 
fails for the general mixed flow problem. We also give many new estimates of the exponent 
9 associated with the exponential tail of the angular tumbling time distribution in the case of 
simple shear. 



1. Introduction 

Advancement in the techniques of fluorescence microscopy have led to direct observation and 
experimentation of single polymer chains (like a A- DNA molecule) in fluid flows in the recent 
years [Il[2l[3llll[5l[6l[71[8]. Several works on theory and simulations have also addressed this 
problem, mostly for the case of "simple shear" flow [9l [lOl [TTl [121 1131 ES ES]- This problem 
has biological relevance jl6j and in particular plays a crucial role in the blood clotting process 
|17j . But planar fluid flows need not be simple shear only, and more complicated flows which 
are either predominantly elongational or rotational have also been studied O ITHJ [191 [20] . 

In general, any planar flow of the form v = VxX + Vyij can be considered as a linear 
superposition of a rotational component with vorticity oj = — and an elongational 

component with strain rate e = + [H]- The "mixed flow" refers to an intermediate 

flow between pure elongational flow {oj = 0) which creates large chain deformation, and 
pure rotational flow (e = 0) which is responsible for angular rotation. "Simple shear flow" 
is the limit where the magnitude of rotational component exactly equals that of elongatinal 
component i.e |a;| = |e|. In this type of flow a polymer tumbles as well as undergoes coil-stretch 
deformation. In Cartesian coordinates the force vector in a mixed flow can be represented as: 
f(t) = (7y(i), a72;(t), 0) where 7 is a constant rate. The parameter a G [—1,1] controls the 
relative magnitude of elongational and rotational components. In the limit a — )• this flow 
reduces to simple shear flow in the X direction. On the other hand, the two extreme limits 
a — )• -|-1 and a — )■ — 1 correspond to pure elongational flow and pure rotational flow respectively 
(see Fig. [I]). 

A polymer in a flow field has been modeled theoretically as Rouse chain |21l |22] with FENE 
end-to-end constraint |19j , as a FENE dumbbell [101 [TT] and as a liner (but non- Rouse) chain with 
global contour length constraint [13]. Very recently we have shown |15j that analytically exact 
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Figure 1. Schematic diagrams of different regimes of planar mixed flow, (a) The pure elongational flow 
(q = +!)■ (b) Simple shear (a = 0). (c) The pure rotational flow {a — —1). Solid arrows represent the force 
vector field and dotted arrows represent the flow directions. 



static distribution functions for the polar coordinates of the polymer's end-to-end vector R(t) 
are possible to be derived for a Rouse chain (without any nonlinear constraint) in simple shear 
(a = 0). Furthermore for the latter flexible polymer, for the first time, we presented an analytical 
calculation for the angular tumbling time distribution function. This is very nontrivial as R(t) 
(although Gaussian) is a non-Markovian process, for which in general calculation of first-passage 
properties like tumbling are very difficult. It was shown that the probability density function 
(PDF) of angular tumbling time r, goes as exp(— ^t/tq) (where tq is the polymer's longest 
relaxation time in the absence of any external flow) and 9 in general depends on the shear rate. 
In the limit of strong shear the exponent 6 goes to a constant value and it was calculated (in Ref. 
|15] ) by a method known as "independent interval approximation" (IIA) \23 \ \2A \ \25 \ I26j. The 
method of IIA was applicable in the limit of strong shear as the X-component of R(t) became 
a smooth process (defined later) in that limit. 

In this paper we study the problem of a linear Rouse polymer in general mixed flow. We 
have derived the exact autocorrelation functions for the Cartesian components of R(t). We show 
that the X-component of R(i) is a non-smooth process and hence analytically IIA method is not 
applicable to find the exponent 6 in the general mixed flow. In the special case of "simple shear" , 
we have also investigated this exponent in more detail (than earlier [T2]) and got few comparative 
estimates of constant 9 using various available IIA schemes. For the general mixed flow we have 
numerical estimates of 9 from observed exponential tail in the tumbling time distribution in our 
simulation. 



2. General correlator for polymer in mixed flow 

We study the linear Rouse model |2H [22] of a polymer chain of N beads connected by harmonic 
springs (each has spring constant k) in a planar mixed flow. The position vector of nth bead 
(n = 1, 2, N) at time t is denoted by r„(t). For internal beads (n = 2, 3, A^ — 1) the position 
vector evolves via the Langevin equation : 

dr (t) 

= -k{2vn - r„+i - r„_i) + f„(t) + 77(n, t) (1) 

where the vector fn(t) = {•^yn{t),ct'^Xn(t),QY' denotes the planar mixed flow force field with 
the constant rate 7. The vector r]{n,t) = {rii{n,t),r]2{n,t),r]s{n,t))'^ is the thermal white noise 
with zero mean and a correlator {T]i{n,t)r]j{n' ,t')) = CSij^n,n'S{t — t'), where i,j = 1,2,3 and 
n = 1,2,..., N. The noise strength ( is proportional to the temperature and in the Eq. ([1]) 
the force strengths are scaled by viscosity. It is convenient to introduce dimensionless so-called 



Weissenberg number Wi = 7T0, where tq = A^^/fcvr^ is the longest relaxation time of the polymer 
in the absence of any external flow. Note that the two end-beads (for n = 1 and n = N) feel 
only one-sided interaction of the harmonic spring. The position vectors of these two beads can 
also be included into Eq. (1) if we add two hypothetical beads at n = and n = N -\- 1 with 
the "free boundary conditions": ro(t) = ri(t) and rAr+i(t) = r]\f(t). 

In the large limit, the discrete index n can be treated as a continuous variable u and Eq. 
([!]) then becomes 

dr(u,t) , d'^riuA) ^, , . 

— ^ = k ' + i{u, t) + r?(n, t) (2) 

with the free boundary conditions : ^'"g"'*^ U=o = and \ u=N = 0. In this continuum limit 

the noise correlator and the force field are given by t)r]j{u' , t')) = (6ij6{u — u')5{t — t') and 

f(ii,t) = (7?/(w, i), a7x(u, t), 0)'^ respectively. 

We solve Eq. ([2]) by Fourier cosine decomposition using the transformation : 

1 TT 

fp(t) = — / COS (^^^ r{u,t) du p = 0,1,2,... (3) 

where each mode evolves independently. For the mode coordinates the problem reduces to the 
single particle problem (i.e a single bead connected to the origin by a harmonic spring in a 
planar mixed flow) with the modified spring constant kp = /tq and noise strength C,' = C,/2N: 

^lM = _kpVp{t) + Ut) + f,p{t) (4) 

where fp(t) = (7yp(t), a75;p(t), 0)^ and r)p(t) are transformed force vector and white noise 
respectively. Now Eq. (jj]) which are coupled equations of the mode coordinates yp{t) 
and Zp{t) can be solved using the method of Laplace transform and we obtain the exact two- 
time correlators of the mode components. 

Our main quantity of interest is the end-to-end vector defined as: R(t) = r{N,t) — r(0, t). 
Using the inverse Fourier cosine transform we find that the correlators associated with the 
Cartesian components of R(t) are related to the correlators associated with the Cartesian mode 
components as below (for finite time increment r > 0) : 

00 

{Ri{t)Rj{t + T)) = 16 i^ipit)rjpit + t)) (5) 

p=l,3,5,... 

where Ri and fjp are the Cartesian components of the end-to-end vector and the pth mode 
vector respectively. From the above equation various correlation functions associated with the 
Cartesian components of R(t) can be obtained explicitly. 



2.1. Exact autocorrelation functions of the end-to-end vector 

Using Eq. dSJ we have derived the exact autocorrelation functions of the Cartesian components 
of R(t) in the stationary state limit (t — t- 00) with a finite time increment r > 0, both in the 
a > and a < regimes separately. 

Elongational regime (a > 0): In this case the elongational effect dominates over the rotational 
effect. We note an important point that in the elongational regime the asymptotic correlators of 
R(t) exist only under the condition: > Wi-^a (for p = 1,3, 5, ...). Thus the correlators diverge 
for the lower modes and hence one would expect huge instability in the polymer's configuration; 



it will get infinitely stretched and will not attain a stable state. We give the results in the limit 
t — 7- oo below, and express them as functions of a scaled time-increment variable T = t /tq : 
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Rotational regime (a < 0): In this regime the amount of vorticity is more dominant than the 
elongational component. We find the correlators in the limit t — )• oo to be the following : 
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Note that these correlators (Eqs. [T2]-[I7]) are in fact the analytic continuation of the results in 
Eqs. [6] -[11] for a > to a < 0. We also note that from the above list of the correlation functions 
we can obtain any static correlator by setting T = in relevant equation. To investigate any 
dynamic properties (like tumbling) we need to keep T 7^ 0. 



2.2. Static distributions 

Since the Cartesian components of R(t) = {Rx{t), Ry{t), Rz{t)Y' are all Gaussian variables, their 
joint PDF is also Gaussian and given by: V{Rx,Ry,Rz) = (27r)-3/2|c|-i/2 exp(-iR^CR), 



where C denotes the covariance matrix involving the static correlators. We also can 
obtain the joint PDF of the polar components of R(t) using the standard transformation: 
V{R,'d,(l)) = P(i2cosi?coS(^, i?cos?9sin(/), iZsint?) x iJ^cosi?. In particular the azimuthal angle 
((/)) distribution Q{(j)) provides some physical insight. Integrating V{R,'&,4>) over R we get the 
joint angular PDF S{^, (j)) and again integrating it over i? we finally get the ^-distribution, 
In the elongational regime (a > 0): The covariance matrix is: 




(18) 

and the determinant is |C| = a/, where / = ad — b'^; for the exact expressions of d, ...,c see 
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Appendix A Finally the azimuthal angle distribution is 

|C| 1 



27rc-v/7 (a cos^ (p + d sin^ <p — b sin 20) 



(19) 



The most probable value of azimuthal angle from the above equation is given by : (pm = 

itan-i (^)- Now, in the pure elongational limit a — t- +1 we get d ^ a (see [Appendix A[ ) 

and therefore (pm — ?• '7r/4. So the polymer stays along the </) = 7r/4 plane for most of the time 
in this case which is expected intuitively. We can calculate the full width at half maximum 
A(j) of Q{(j)) near (pm ~ '/r/4 and it is given by : cos(2A(^) = 2 — (a/6). For Wi — oo, we get 
a/6 ~ 1 and hence A(p — )■ 0: so the Q{<p) is sharply localised which is quite expected since the 
pure elongational flow is vorticity free where the polymer hardly tumbles. But such a limit is 
not actually realisable as Wi — )• oo violates the mode stability condition > Wi^/a, and so 
such a aligned polymer will experience an immediate catastrophic stretching. 

In the rotational regime (a < 0): In this case we write the covariance matrix as : 

(20) 




and again the explicit values of a, b, c and d are given in [Appendix A The azimuthal distribution 



is : Q{(p) = ( — ^^7^ ) /(acos^ (p + dsin^ (p — bsm2(p). In pure rotational limit a — )• —1, we get 

\2^TC^/f J 

dK, d and 6—7-0, hence we have Q{(pi) = l/27r =constant. Clearly, the pure vorticity dominated 
limit is quite expected to have a flat (p distribution. Moreover as a — )■ —1, dK, a k. c and 6 — )• 0, 
leading to all eigenvalues of the matrix C becoming identical. This means that the polymer 
approaches a spherical shape as a — )■ — 1 (pure rotational flow) . 



2.3. Dynamic properties 

Now we concentrate on the first-passage question of the "tumbling" process of a polymer in 
mixed flow. The tumbling event is either defined as the successive occurrence of coiled states 
(radial tumbling) [3 HHj or as successive crossing of a fixed plane, say = 0, by the vector R(i) 
(angular tumbling) [5l[6l[T0]. The definition of radial tumbling is subjective, since a tumbling 
event is marked by |R(i)| falling below a value Rq that defines the coiled state. But the angular 
tumbling time r is objective and is simply the interval between two successive zero crossings of 
the process Rx{t)- The PDF of r asymptotically goes as P{T) ~ exp(— 6^), where T = t/tq is 
the scaled time- increment. 

The mean density of zero crossings for a Gaussian stationary process is given by: p = 1/(t) = 
y^—A"{0)/ir, where A{t) is the normalised correlator of that process and r denotes the intervals 



between zeros. The quantity p is finite when ^'(0) = and j4"(0) is finite. Such processes 
are called "smooth" and for them one can apply IIA to find the nontrivial exponent of the 
asymptotic distribution of the intervals between zeros |23t [M] . 

In the present case we are only interested in angular tumbling of the polymer. In the 
elongational regime (a > 0) angular tumbling is very rare event in the limit a — )• +1 due to 
sharp localisation of (/(-distribution for which (p hardly changes sign. To investigate the angular 
tumbling for a < we proceed with the autocorrelation function Cxx{T) = \\m.t^^{Rx{t)Rx{t + 
tqT)) (see subsection 2.1) and consider the normalised correlator A(T) = Cxx{T)/Cxx{^)- To 
know the small T behavior of A{T) we first obtain the exact Laplace transform of CxxiT) : 
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where /3 = Wiy/a. From the above equation, the expression of A(s) in large s limit (i.e s — )■ oo) 
can be obtained easily and inverting that we derive the expression of A(T) for T — )• as below: 
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where Cxx{0) = {Rx) ^^'^ static correlator given in Appendix A It is clear from the above 
equation that both ^'(0) and ^"(0) diverge. Thus Rx{t) is a "non-smooth" process in the sense 
that the density of zero crossings is infinite and hence the method of IIA can not be applied. 

Although IIA is inapplicable, we have simulated the Langevin dynamics of Eq. ([1]) and 
curves for P{T) were obtained (see Fig. [2]). The exponential tails fall more sharply {6 increases) 
with increasing negative a. Furthermore the appearance of a peak (for a < 0), signifying a 
characteristic cyclic tumbling, gets more prominent with increasing negative a. 



3. Special case of simple shear flow 

The limiting case of "simple shear" (a = 0) is the only solvable case for P{T) using the method 
of IIA [15], that too for large shear rate Wi — )• oo. In the limit a — )• 0, Eq. (|22p yields 
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Now we note that in the limit of strong shear (Wi — )■ oo) the singular contribution of the second 
term proportional to \/T in Eq. ()23|) is ignorable. In the limit of T — >• 0, Eq. ()23p has the 



limiting form A{T) ~ 1 — (GO/vr )r as T — )• 0. Hence, in the limit Wi — )• oo the process 

Rxit) becomes "smooth" [25] and IIA can be applied to find the constant 9. Clearly we have 
^'(0) = and ^"(0) = — 120/vr. Thus we get finite density of zero crossings p = \/120/vr^. If 
we assume P{T) to be exponential not only asymptotically but also over the full range of T, the 
numerical value of p gives a crude estimation of 9 as below: 

9^ p 0.353298 •• • 

We note that the expression of A{T) for any T in the limit of Wi — )■ oo is : 




Figure 2. (color online) Linear-log plot of P{T) versus T. Values of the parameter a are respectively: a — 
(the upper curve, with red dots), a = —0.1 (the middle curve, with green dots) and a — —0.3 (the lower curve, 
with blue dots). The fitted exponential function for the simple shear case (q = 0) is : 0.00012 x exp(— 0.324r) 
(where 9 is used from method IV in section 3). All data are for A'' = 10 and 7 = fc = 1. 
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Now we show below different estimates of 6 from various IIA calculations. 

(I) First we consider the "clipped" variable xi^) = sgn{Rx). Let r{T) = ^^^^^^(^^2^'^^'^ be the 
normalised asymptotic {t — )■ 00) correlator of xii)- One can express: r(T) = X]^o(~-'-)"P"(^)' 
where Pn(T) is the probability that the interval T has n zeros. The strategy is to write the 
joint PDF of n successive intervals as the product of the distribution of single intervals by IIA 
|23^ I26j. Then one ultimately finds a relation between the Laplace transforms of P{T) and r{T) 
(which are denoted by P{s) and f(s) respectively) as below: 



Pis) = [1 - (p/2)s(l - sfis))]/[l + ip/2)sil - sfis))] 



(25) 



where p is the density of zero crossings. The asymptotic PDF P{T) ~ exp{—9T) shows that 
P{s) has a simple pole at s = —6 i.e the denominator in Eq. (|25|) has a simple zero at that point. 
Now, the correlator of the clipped process is given by: r{T) = (2/7r)sin-i(yl(T)) (this relation 
holds for any Gaussian process [27]), where A{T) is given by Eq. (|24|) . Using this relation in 
Eq. (|25|) and equating the denominator to zero we get (with a transformed variable x = e~'^): 
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Solving the above equation numerically for s = —8 we obtain : 

6* « 0.356203 •• • 



(26) 



(II) Two other IIA schemes are given in Ref. [23]. Let U (T)dT be the conditional probabihty 
that a zero of occurs between (t + T) and {t + T + dT), given that there is a zero crossing of 
x{t) at time t. It has been shown by J. A. McFadden that P{s) and U (s), the Laplace transforms 
of P(T) andC/(T) respectively, are related by : P{s) = U{s)/[1 + U{s)]. The expression of C/(T) 
for a Gaussian process is given by [23] : 

where 

r" = (2/^)M23(l - A2)-3/2 . H = M23 {Ml - M^y^'^ 

M22 = -^"(0) (1 - yl^) - ; = A" {I - A^) + AA'^ 

The dependence on T is understood in the above equations. We note that the function U{s) 
raises some computational difficulty as it does not exist at s = 0, because U (T) does not vanish 
as T — >• 00. Also one can observe that as T — )• 00, U{T)dT — )• pdT. If a new variable is defined 
as [23]: = U{T) — p, with Laplace transform V{s) = C/(s) — (p/s); then U{s) can be 

substituted by V{s) and the computation can be performed easily. Setting the denominator of 
P{s) to zero we get the equation : s + p + sV{s) = 0. Solving this equation numerically for 
s = —6 we have : 

e ^ 0.257136- •• 

(III) In the third scheme, instead of assuming that all successive zero-crossing intervals are 
independent, it was assumed in Ref. [23] that a given interval is independent of the sum of 
the (2m + 2) earlier intervals (for m = 0,1,2,...). This is called "quasi independent interval 
approximation" (QUA). By QUA one gets : P{s) = {U{s) + g{s))/[2 + U{s) - g{s)], where 
g{s) = L{r"{T)/4p}. Again solving numerically for the zero of the denominator of P{s) (using 
the newly defined variable V{s) instead of U{s)) at s = —9 we get : 

e ^ 0.213558- •• 

(IV) In Ref. [15] another approximation r{T) w A(T) was used that gives the value of 6 
very close to the numerics (see a = case in Fig. [2]). This approximation yields : 

6* « 0.323558 •• • 

We see that various estimates of 9 differ slightly, but upto first decimal place (method III 
being an exception) the number is roughly ~ 0.3 . 

4. Conclusion 

We have exactly calculated the most general spatio-temporal correlation functions for a Rouse 
Polymer in planar mixed flow, in the steady state. Using the exact correlators of the mode 
coordinates, we explicitly write down the various correlators associated with the Cartesian 
components of the end-to-end vector R(t). The static correlators may be used to find the 
static probability distributions associated with the polar (??) and azimuthal {(p) angles of R. 
These static correlations give important information about the average shape of the polymer. 

The more important concern in our paper has been the dynamics of angular tumbling of the 
vector R(t) in mixed flow, and the distribution of time scales associated with these flrst passage 
events. We find numerically that for rotationally dominated flows (a < 0), the distribution has 
a peak and then an exponential tail for large times. We demonstrated that analytical estimates 



of the decay constant of the exponential tail can not be estimated using the method of IIA in 
Q < cases, as the process Rx{T) is not smooth for any limiting value of Wi. 

The method of IIA was successfully applied by us [15] in the special case of "simple shear" 
(a = 0), using a particular type of approximation (namely method IV in section 3). Here 
we have given many new analytical estimates of the number 9, using many different types of 
approximations following the classic work of Ref. [23]. Although the numbers are slightly 
different from each other, up to the first decimal place 9 « 0.3. 
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Appendix A. Static correlators 

In this appendix, we give the exact expressions of the static correlators, which we get by putting 
T = in the equations ([6|) to (fH 



The static correlators for q > 
Using the standard series sums 
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to get : 
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Similar steps yield the following : 
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where x = y/wla^^^. 

The static correlators for a < : 
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c- = {Rl) = (A.10) 

where y = ^JWl/2\a\^^^ . We checked that in the Umit a — )• 0+ as weh as a — )• 0— the above 
static correlators smoothly go to the exact expressions of the static correlators in "simple shear" 
that are published in Ref. [T5] . 
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